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ABSTRACT 



The dynamical evolution of the globular cluster systems in galaxies is predicted, 
based on the standard dynamical theory normalized to the example of the Milky Way. 
The major processes varying with the galactocentric distance are the tidal shocks and 
dynamical friction. Our simple model explains, on a quantitative basis, the observed 
differences of the inner and outer populations of globular clusters. We can thus calculate 
corrections for dynamical evolution for the luminosity function of globular clusters with 
the assumption that the initial luminosity function is identical in all galaxies (and we 
can test this assumption as well, in certain cases). Then we can compute the expected 
distribution of absolute magnitudes and compare it with the observed distribution of 
apparent magnitudes to estimate the distance moduli for M31 and M87. Using this 
new method we find (dmu3i = 24.05 ± 0.23, dmyisi = 30.83 ± 0.17) as compared 
to current best estimates using other methods of (dm,M3i = 24.30 ± 0.20, dniM87 = 
31.0±0.1). As a check on the method we compute, and compare with observations, the 
differences between the inner and outer globular clusters in all three galaxies. This new 
method, coupled with HST observations, promises to provide an independent method 
of estimating distances to galaxies with recession velocities ^ 10,000 km s , or d 
lOO/i" 1 Mpc. 



Subject headings: globular clusters: general — galaxies: star clusters — galaxies: dis- 
tances and redshifts — galaxies: individual (M31, M87) 



1. Introduction 

Evolution of the globular cluster systems (GCS) in our Galaxy and external galaxies has 
received a burst of recent interest. New HST observations of the GCS in the Virgo ellipticals allow 
a direct comparison of the globular cluster luminosity functions (GCLF) in different galaxies. The 
assumption that the peak of the GCLF occurs at a constant luminosity for galaxies of a given type, 
and varies slowly and in a predictable manner with the galactic type, has been used widely as a 
distance indicator ( [Harris 1991| ; Jacoby et al. 1992; Harris 1996| ) although the theoretical basis for 
this expectation remains unclear. 
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On the theoretical front, Gnedin fc Ostriker (1997) and Murali fc Weinberg] (1996a, b,c) extend- 
ing earlier work of Aguilar, Hut, & Ostriker (1988) and Chernoff &: Weinberg (1990)| have shown 
that significant evolution of the GCS is to be expected over the Hubble time. The inner popula- 
tions of globular clusters are most affected by the tidal shocks and dynamical friction. The first of 
these effects discriminates preferentially against the less dense and less massive clusters, the latter 
affecting the most massive clusters, to shape the currently observed distribution of luminosities. 
We have shown (Gnedin 1997) that the expected statistical "brightening" of the inner clusters (due, 
really, to the destruction of some of the less luminous clusters) is indeed observed in the three best 
studied galaxies, the Milky Way, M31, and M87 (see also [van den Bergh 1996 ). In this paper we 
examine the hypothesis that a common initial distribution of clusters existed in the inner and outer 
parts of all the studied galaxies with the apparent variations due to the influence of dynamical 
evolution. And then we show how the method can be used to obtain an independent estimate of 
the distance for galaxies with an observed distribution of globular cluster apparent magnitudes. 



2. The mass — density distribution in the Milky Way 

The primary processes driving globular cluster evolution are two-body relaxation, evaporation 
of stars through the tidal limit, tidal shocks, and dynamical friction. Relaxation is the only internal 
effect and in zeroth approximation does not depend on the position of the cluster in the galaxy. The 
tidal effects do depend on the density of the galactic environment and tend to destroy less dense 
clusters. The rate of inspiral due to dynamical friction is proportional to the cluster mass, M, thus 
preferentially destroying the most massive clusters. Both the tidal effects and dynamical friction 
are strongly enhanced in the inner part of the galaxy and are weak on the outskirts. Therefore 
it is possible that the outer populations of globular clusters may still bear the form of the initial 
distribution, whereas the inner population is affected by the extrinsic dynamical effects. Passive 
evolution, through stellar interior processes and two-body relaxation, is expected for all clusters. 
A more careful treatment would allow for the interaction between internal and external effects. 

First we shall make an attempt to reconstruct the original distribution using the outer cluster 
sample. The simplest generalization of the GCLF is a two-dimensional distribution of the cluster 
mass and density. A useful estimate of the cluster density is p = M/Rl, where R h is the half- 
mass radius. Since only the Galactic clusters are well resolved, we use the sample of the Milky 



Way clusters tabulated by Djorgovski (1993) . The core-collapsed clusters were eliminated from 
the sample due to the observational difficulty of estimating their radii. We also rejected the few 
very distant clusters (beyond 60 kpc from the Galactic center) as being of uncertain origin. The 
remaining sample consists of 94 clusters. We set the boundary between the two samples at the 
galactocentric radius R\/2 = 5 kpc, thus putting 36 clusters in the inner and 58 clusters in the 
outer groups. The outer sample has more members because we expect the inner population to 
be significantly depleted. This choice also maximizes the difference between the two populations 
within the range 4 < R\/2 < 7 kpc. 
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Fig. 1. — Distribution of the inner (filled dots) and outer (stars) Galactic globular clusters on the 
log M — log p plane. The solid lines show the adopted initial distribution function (eq. [H) where the 
ellipses are la, 2a, and 3a levels, and the small circle is at the center of the distribution. The dotted 
lines mark the region allowed for the inner clusters by dynamical processes: above the horizontal 
line dynamical friction is important, to the left of the vertical line tidal shocks are important, and 
below the diagonal line relaxation will lead to core collapse and subsequent disintegration of the 
clusters. The equivalent of the first two of the survival boundaries for the outer clusters are not 
shown as they do not significantly constrain cluster properties. 
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Figure |l] shows the distribution of the defined sample of Galactic globular clusters in the 
logM — log/9 plane. Throughout, we use an assumed mass-to- visual-light ratio of MjLy = 3 to 
estimate masses when velocity dispersions are not available. The area occupied by the inner clusters 
is visibly constrained by the survival boundaries of dynamical friction and the tidal shocks. The 
boundaries for the inner clusters (the horizontal and vertical dotted lines) are defined as the median 
clusters that would marginally survive in the age of the Galaxy conservatively estimated as to = 12 



Gyr. For details see |Gnedin fc Ostriker (1997) 



As noted by |Aguilar, Hut, fc Ostriker (1988) and |Gnedin fc Ostriker (1997) , dynamical friction 



does not pose a serious cutoff on the distribution of the clusters in Milky Way. It will, however, 
be seen to be important for the M87 globulars because their distribution function is much better 
sampled (see §||). The tidal shocks are, on the contrary, very important in destroying the low 
density clusters. 

Notice a very strong correlation between the mass and density of the outer clusters. The 
Pearson correlation coefficient is 0.76 and the Spearman coefficient is 0.77, so the probability of no 
correlation is less than 10 -16 . The two-dimensional least-squares fit gives a useful relation between 
the half-mass radius and mass of globular clusters that has not been investigated before: 

/2xlO 6 M \°- 63 

R h = [ M ° j PC (i) 

or approximately Rh oc M -2 / 3 . It is worth noting that this relation predicts a very steep dependence 
on mass of the tidal shock destruction time, t s h oc MRT 3 oc M 3 , and a shallow dependence on mass 
of the relaxation time, t r h oc Af 1 / 2 /?^ 2 oc Af" 1 / 2 . As to whether or not this correlation is primordial 
remains to be seen, since the present sample could have been affected by relaxation itself. 

Including the correlation in the bivariate distribution reduces the occupied space on the log Af — 
logp plane and maximizes the distribution function /(log Af, log p). While the log-luminosity dis- 
tribution is known to be approximately Gaussian, the distribution of log p is certainly skewed. 
Nevertheless, we assume the bivariate Gaussian shape for the combined distribution for the sake 
of tractability. Thus the adopted primordial distribution function of globular cluster properties is 
taken to be 



f(x, y) dxdy = ° + ^ exp 

ATXO\0~2 



{y — ax — b) 2 (y + x/a — cY 



2aj 2a\ 



dxdy, (2) 



where x = log/3, y = logM; a = 0.35 is the slope of the correlation, and the constants b and c 
are fixed by the center of the two-dimensional distribution. The dispersion along the line of the 
correlation, ai = 3.70, is much larger than that in the perpendicular direction, a\ = 0.32. The 
units adopted are such that the cluster mass is in solar masses and the density is in M Q pc" 3 . The 
solid lines in Figure [l] show the adopted distribution function. It provides a satisfactory fit to the 
outer clusters (denoted by stars). 
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Given the initial distribution, we can evolve it in time and compare with the present data for 
the inner population. The most straightforward approach is to assume a hard cutoff at the density 
p s h and mass M^f such that the clusters less dense or more massive than that would have been 
destroyed within the Hubble time by the tidal shocks or dynamical friction, respectively. Such 
survival boundaries are shown on Figure |l]. Then, since for most observed GCSs we do not know 
individual values of and hence p, we integrate the evolved distribution function over p (hence x) 
to obtain the evolved mass distribution which can be used to compute the luminosity distributions 
when studying the GCSs in external galaxies. The result of the integration is 

F(y) dy = — L_ e ~(^o W 0(p - y) Jerfc(^) dy, (3) 

V I'KOy * 

, ,— y — c a(y — b) a 2 1 

w(y) = ^ag+- — - 2 , a = — ^ + 2 , 

where p = logM^, g = log p s h, 9(p — y) is a step function, and erfc(tn) is the complementary 
error function. Equation (^) should also be renormalized to the size of the remaining sample. An 
important consequence of the modified distribution is the shift of the peak, given by the nonlinear 
equation 

Ay = y - y - - ( — ) a 2 ^— (4) 

\ dy J y y/n erfc(-07) ' 

where we neglected the dynamical friction correction as it is usually not important for the bulk of 
the sample. This equation can be easily solved iteratively to obtain Ay. When all parameters of 
the initial distribution (eq. ||) are fixed, the predicted shift is a function only of g, the amplitude 
of the tidal shocks. It was calculated approximately, equating the inverse destruction rate of each 
cluster due to the shocks, t s h (see Gnedin & Ostriker 1997), to the Hubble time. A correction factor 
of order unity can be used to calibrate g against the data. Taking t s h/0.3 as true destruction time 
gives the observed shift between the inner and outer clusters in Milky Way (Ay ~ 0.16), and we 
adopt this factor here. The correction in magnitudes is related to the shift Ay by Athq = — 2.5Ay. 
The dispersion of the inner sample can also be calculated similarly. 

The observed and predicted mass distributions of the Galactic globular clusters are shown 
on Figure ||. The agreement is excellent. Our correction reproduces well both the mean and 
the dispersion of the inner sample. The first of these two checks is of course not significant (but 
only a check of consistency) since we fixed the factor 0.3 by the requirement that the observed 
and computed inner peaks agree at My in = —7.73. But the correct predicted reduction of the 
dispersion Act = 0.25 ± 0.12 mag seen between the two samples (Aa i, s = 0.23 ± 0.16) provides a 
fair test of the hypothesis that the inner population had the same initial properties as the outer 
population but was modified by dynamical effects. Table [l] provides a comparison of the predicted 
and observed changes in the mean and the dispersion, as well as the associated uncertainty. 
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Fig. 2. — The histograms of the inner (solid line) and outer (shaded) samples of the Galactic 
globular clusters. The thick dashed line and the thin solid line are the Gaussian fits to the inner 
and outer populations, respectively. The thick solid line is the predicted evolved distribution (eq. 
|J) for the inner clusters assuming that they started with the same properties as the outer clusters 
but suffered from dynamical evolution. 
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3. Predicted evolution in M31 and M87 

The predicted evolution of the mass (and luminosity) function of globular clusters can now be 
immediately applied to the samples in M31 and M87 to test if the predicted values of both Amo 
and Act match observations for these independent samples. We use the same data as Gnedin (1997) 
who gives the appropriate references. 

The amplitude of the tidal shocks scales with the density of the galaxy at a given position, 
which is p oc M ga i(R) / R? oc v 2 /R 2 for a spherically symmetric system. Here v c is the rotation 
speed (or the velocity dispersion for elliptical galaxies) at the median radius R me d of the globular 
clusters from the center of the galaxy. While the velocity can be directly measured, only angular 
distances are available for external galaxies so the physical distance will scale with the assumed 
distance to the galaxy, R me d oc D. The dynamical friction cutoff also depends on both v c and R„ ie d 
through Mdf oc VcR^^, where we neglect the slowly varying logarithmic factor. Even though both 
p s h and Mdf depend on the distance D to the galaxy, the correction (eq. 0]) can be applied in a 
self-consistent manner without any initial assumptions about the distance. The method proposed 
below is fully independent and will predict, rather than use, the distance D. 

If the observed samples is rich enough it could be divided into inner and outer halves, or 
even further; the method can be applied to each of the subpopulations separately and the distance 
estimated as average over the individual results. In principle, the narrower the division of the sample 
into radial bins the better the accuracy of the median estimator R me d', in reality, the subsamples 
should have a large enough number of members to represent well the luminosity function. For 
M31 and M87 we take the inner and outer halves of the sample, as defined by Gnedin (1997). For 
M31, the inner clusters have the median radius at Ri n ,med = 13/5, whereas the outer clusters have 
Rout,med = 50'. The rotation velocity of M31 is estimated to be v c ~ 300 km s _1 for 10' < R < 20', 
and v c ss 330 km s^ 1 for R > 20' QFederici, Fusi Pecci fe Marano 1990[ ). 

Let us demostrate how the method works on the example on the inner population of globular 
clusters in the Andromeda. The mean magnitude of the subsample mo,i n = 15.93 (Gnedin 1997) 
so the trial value of the distance modulus is taken to be dm^ = niQ^ n — My = 23.26, where 
My = —7.33 is the center of the assumed initial distribution of globular clusters. Using this 
distance we calculate the critical shock density p s h,M3i = 5.4 x p s h,MW an d the corresponding 
correction to the peak magnitude, Amo in = 0.80. The peak of the inner GCLF therefore becomes 
dimmer and the second estimate of the distance modulus is dm^) = 24.06. Now second and further 
iterations are performed to calculate new p s h and to obtain new estimate of dm until convergence 
is achieved. Usually 4 or 5 iterations are sufficient to obtain the distance modulus to an accuracy 
better than 0.01 mag. In case of the inner clusters of M31, the final predicted shift of the peak 
magnitude is Amo,^ = 0.65 and the distance is dm,M3i = 23.91. The same operations are then 
repeated for the outer clusters as well, starting with their mean mo >ou t = 16.71. The results are 
given in Table |]. Now combining the two predicted changes for the inner and outer populations we 
get an estimate of Amo an d Act for that galaxy. The predicted change agrees with the observed 
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values within the errors (Table [l]). 

The distance modulus of M31 has been obtained by a number of other methods. The Cepheid 
variables give dm = 24.43 ± 0.16 flFreedman k. Madore 1990) ), RR Lyrae give dm = 24.34 db 0.15 



( Pritchet fc van den Bergh 1987] ), fitting color-magnitude diagrams of population II giants gives 



dm = 24.40 ± 0.25 flMould fc Kristian 1986| ), and novae give dm = 24.04 ± 0.20 ( |Cohen 1985| ). The 



average of all these methods indicate dmyizi = 24.30 ±0.20. Our method predicts dmyizi = 24.05 ± 
0.23 (averaged over the inner and outer samples; see Table |2|), in good agreement with the other 
methods. Note that our result was derived independently from any other distance determinations 
to the galaxy. 

In case of the Virgo giant M87 the dispersion of radial velocities v c is used instead of the 



rotation velocity as an estimate of the galaxy mass. Merritt fc Tremblay (1993) give the mass 
of M87 of 6 x lO 12 Af at R = 50 kpc, so that v c ss 730 km s _1 . The globular clusters are 
distributed over a larger area; the median radius of the inner population is at Ri n ,med = 2. '4, and 
of the outer population is at R ou t,med = 5/1. Repeating the same steps as for the M31 sample, we 
obtain a very good agreement of the predicted changes with the observed Amo = 0.26 ± 0.03 and 
Act = 0.04 ± 0.03 (see Table [l]). The distance modulus to M87 is essentially the same as inferred 
from the inner or outer samples and is dmyi87 = 30.83 ± 0.17, in remarkably good agreement with 
the value of dmyi&7 = 31.0 ± 0.1 as obtained for the center of the Virgo cluster based on Cepheids, 



planetary nebulae luminosity function, and surface brightness fluctuations (quoted by Harris 1996 ). 
The largest source of error in the method proposed here is due to the uncertainty in the peak of 
the Galactic globular clusters, 8 My = 0.16 mag. 

Because of the extreme richness of the clusters in M87, the effects of dynamical friction should 
be detectable. Indeed, if we keep the same division of the clusters into the inner and outer pop- 
ulations but consider only the brightest tails of both, with V < 20.9, the outer clusters become 
relatively brighter (as the brightest inner clusters are depleted). This effect is illustrated in Table 
||. Statistically significant, it shows lack of very bright, and presumably massive, clusters in the 
inner region, as is expected, if they spiraled into the center by dynamical friction. The medians 
provide an even stronger case at the 5a level. 

Clusters at R = 12 kpc would reach the center of M87 in the Hubble time due to dynamical 
friction if they are more massive than 1.4 x 10 7 Mq flBinney fc Tremaine 19871 , P- 428). This 



corresponds to the visual magnitude of Vjf = 19.3. In accord, only 4 clusters brighter than Vdf 
are present in the inner sample while the outer sample has 9 as bright members. As suggested by 



Ostriker, Binney Saha (1989) , Pesce, Capuzzo-Dolcetta, fe Vietri (1992) , and Capuzzo-Dolcetta 



(1993) , the dynamical friction rate is strongly enhanced in triaxial halos with box orbits; this effect 
can lower even further the mass cutoff. Therefore dynamical friction is important in determining 
distribution of very bright clusters in M87. 
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4. Discussion 

By applying the corrections for dynamical evolution, the GCLF can be used as an improved 
distance indicator. For example, the method can be used for determining distances to the Virgo 
cluster and, possibly, to the Coma cluster. HST offers the unprecedented opportunity to detect 
faint globular cluster systems around galaxies in the Virgo cluster. Having a large sample of such 
galaxies is essential to eliminate another source of error, the difference in metallicity of different 
galaxies (see Ashman, Conti fc Zepf 1995| ). 



We have hitherto assumed a fixed mass-to-light ratio whereas it is known that both the lumi- 
nosity of individual stars of fixed mass and the distribution of masses within clusters are a function 
of the cluster mean metallicity. An ab initio correction for such effects would be difficult and uncer- 
tain. But an empirical approach which utilizes the fact that for an evolved population, the mean 
color is a good proxy for metallicity, should be possible. The minimization of the effects due to 
metallicity differences could be done as follows. Assume that the (B — V) color of the integrated 
galaxy light is a fair indicator of the metallicity. Then the mean, mo i, of the GCLF of each galaxy 
i is assigned a metallicity correction of the form Srrii = a 5(B — V)i = a[{B — V)i — {B — V)], where 
a is a constant to be determined, in addition to the correction Am^^ for dynamical evolution; 
here (B — V) is the mean color of the Galactic sample. The true mean point is thus 

frii = m ,i + Am dyn> i + a 5(B - V)i. (5) 

We now consider a sample of galaxies all approximately at the same distance (for example, in the 
Virgo cluster) but with different mean colors and metallicities. Then 

(m) = {m ,i + Amdyn,i) +a{5(B- V)»). (6) 
We then minimize the dispersion of the individual means rn~i around the average (m), to determine 



a: 

_ (moj, + Am dyn)i ){5(B - V)j) - {(m 0ti + Am dynti )5(B - V)j) 

(S(B - V)f) - (5(B - V)i) 2 ■ 1 j 

Using this value of a, we obtain a statistical estimator (eq. ||) of the peak of the GCLF to be 
calibrated against the known absolute magnitude, My, of the GCLF in the Milky Way and M31. 
The proposed corrected estimator has the advantage of being justified by simple, relatively well 
understood physics. When used in a real survey of globular clusters in external galaxies, it may 
significantly reduce the scatter among the results for different galaxies. 



5. Conclusions 

The differences in the inner and outer populations of globular clusters in our own Galaxy 
(Gnedin f997) can be explained in terms of simple dynamical processes such as the tidal shocks 
and dynamical friction. Assuming that the primordial distribution of globular clusters is common 



-10- 



for all galaxies, we derive the correction to the turnover point of the GCLF in any galaxy based on 
the dynamical calculations for the Galactic clusters by |Gnedin fc Ostriker (1997)| . The predicted 
shifts between the mean points of the inner and outer cluster populations agree qualitatively with 
the observed differences within the errors. The shifts in the dispersion of luminosity between the 
inner and outer samples of the Galaxy, M31, and M87 are predicted accurately. A slight mismatch 
of the mean points for the Milky Way and M31 of the order 0.2 mag still remains, but is within 
the statistical errors of the methods applied. 

The corrections for dynamical evolution thus allow one to use globular cluster systems in 
external galaxies as an independent distance indicator. For a large sample of galaxies at the same 
distance, as in the Virgo and Coma clusters, the metallicity corrections could also be used with the 
aim of reducing still further the scatter in the derived distances. An extensive survey of the GCS 
in distant galaxies, performed with the HST, may help produce a reliable distance ladder on the 
scale of a hundred Mpc. 



This project was supported in part by NSF grant AST 94-24416. 
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Table 1. Predicted and observed differences in the luminosity function 





Milky Way 


M31 


M87 


Arrio, P redicted (mag) 


0.41 ±0.16 


0.51 ±0.18 


0.25 ±0.17 


Am , o & se ™ed (mag) 


0.40 ± 0.23 


0.78 ±0.12 


0.26 ±0.03 


^predicted (mag) 


0.25 ±0.12 


0.16 ±0.12 


0.07 ±0.12 


^observed (mag) 


0.23 ±0.16 


0.24 ±0.08 


0.04 ±0.03 



Table 2. Distance moduli to M31 and M87 





M31 


M87 


This paper 


24.05 ±0.23 


30.83 ±0.17 


Other methods 


24.30 ±0.20 


31.0 ±0.1 



Table 3. The bright tail of the globular cluster distribution in M87 



Sample N mo (mag) a (mag) fi (mag) 



Inner 108 20.82 ± 0.05 0.63 ± 0.03 20.56 ±0.03 

Outer 93 20.75 ± 0.07 0.85 ± 0.05 20.37 ±0.10 

difference -0.07 ± 0.03 -0.20 ± 0.04 

probability 3.1 x 10~ 2 1.2 x 10~ 5 



Note. - Clusters with V < 20.9 and 1/21 < R < 7'. N is the number of clusters in the 
sample, and mo, c, and [i are the mean, dispersion, and the median, respectively. The data is 
from McLaughlin, Harris, & Hanes (1994). Probability quoted is the probability of Monte-Carlo 
simulations that a given result occurs by chance. 



